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We employ ab initio methods of quantum chemistry to investigate spin-1/2 fermions interacting 
via a two-body contact potential in a one-dimensional harmonic trap. The convergence of the total 
energy with the size of the one-particle basis set is analytically investigated for the two-body problem 
and the same form of the convergence formula is numerically confirmed to be valid for the many- 
body case. Benchmark calculations for two to six fermions with the full configuration interaction 
method equivalent to the exact diagonalization approach, and the coupled cluster method including 
single, double, triple, and quadruple excitations are presented. The convergence of the correlation 
energy with the level of excitations included in the coupled cluster model is analyzed. The range 
of the interaction strength for which single-reference coupled cluster methods work is examined. 
Next, the coupled cluster method restricted to single, double, and noniterative triple excitations, 
CCSD(T), is employed to study a two-component Fermi gas composed of 6 to 80 atoms in a one¬ 
dimensional harmonic trap. The density profiles of trapped atomic clouds are also reported. Finally, 
a comparison with experimental results for few-fermion systems is presented. Upcoming possible 
applications and extensions of the presented approach are discussed. 


I. INTRODUCTION 

Ultracold gases are highly controllable systems ideal for investigating different phenomena of quantum many-body 
physics [1-6] . They can be prepared in a well-defined quantum state, carefully manipulated, and accurately measured, 
and thus can serve as a perfect tool for quantum simulations [7]. On one hand, they can be used to realize highly non¬ 
trivial states of matter. On the other hand, the problems of condensed-matter or other areas of physics can be mapped 
on and solved by such quantum simulators. Reducing the dimensionality of the trapped gases to one dimension brings 
a plethora of new interesting possibilities [8]. Experimental studies of the Tonks-Girardeau gas [9-11] and of the super 
Tonks-Girardeau gas [12] are the first and most eminent examples. 

Levi Tonks (1897-1971) was the first to consider in 1936 [13] the problem of equation of state simultaneously for 
one, two, and three-dimensional gases of hard elastic spheres - this has led him to the concept of a (classical) gas 
of impenetrable particles in ID. Marvin Girardeau (1930-2015) was a real pioneer of the studies of the quantum 
impenetrable gas of bosons, known since then as Tonks-Girardeau gas. In particular, Girardeau investigated in 1960 
the relation between systems of impenetrable bosons and fermions in one dimension [14]. This seminal work attained 
an extreme importance in the cold atom era, and has led to numerous nearly equally seminal generalizations and de¬ 
velopments: studies of ID Coulomb gas [15], studies of trapped bosonic gases in ID [16], studies of quantum dynamics 
and quantum solitons [17], general theory of Fermi-Bose [18, 19] and anyon-fermion mapping [20], investigations of 
super-Tonks-Girardeau gas [21, 22], of ID dipolar gases [23], and much more. From the point of view of the present 
work the most important were more recent generalizations and applications of Girardeau’s approach to soluble mod¬ 
els of strongly interacting ID ultracold gas mixtures [24, 25], and especially to spinor Fermi gases [26-29]. Here we 
investigate precisely the problem of ID fermionic gas of atoms of spin 1/2, trapped in a harmonic potential and in 
the presence of strong interactions. 

Recently, experiments achieving a deterministic preparation of tunable several-fermion systems in a one-dimensional 
trap became possible [30]. This opened the new area of ultracold research on systems with a complete control over 
all degrees of freedom: the particle number, the internal and motional states of the particles, and the strength of the 
interparticle interactions. The fermionization of two distinguishable fermions [31], the formation of a Fermi sea [32], 
pairing in few-fermion systems [33], two fermions in a double well potential [34], and antiferromagnetic Heisenberg 
spin chain of few cold atoms [35] have been investigated experimentally and are a promise of upcoming new and 
equally fascinating research. 
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The aforementioned experiments allowed to observe for the first time the transition between the few-fermion limit 
and the many-fermion limit of trapped atoms at ultralow temperatures. The emergence of the many-body properties 
of the physical systems is crucial across all areas of research. The Fermi gases in one dimension in the many-body 
regime have been studied intensively over the years [36-45] . Recently, several papers used analytic methods and exact 
diagonalization to investigate in detail the few-body regime [46-57]. Nevertheless, new numerical approaches providing 
additional insight on the experimental findings and having the predictive power of proposing new experiments are 
always welcome. Especially methods that can cover both regimes of few and many fermions are of great interest. 

The goal of the present work is to provide new numerical tools to investigate the Fermi systems of few to many cold 
atoms. More specifically, we consider spin-1/2 fermions confined in a one-dimensional harmonic trap and interacting 
via a two-body contact potential. For this aim, we employ a quantum chemistry approach - the coupled cluster (CC) 
method [58-68]. This method has successfully been applied to study various properties of atoms, molecules, and con¬ 
densed phases - see, for instance. Refs. [69, 70] for applications to high-precision spectroscopy of ultracold molecules, 
Ref. [71] for simulations of liquid water properties, or Ref. [72] to determine the structure and characteristics of 
molecular crystals. 

In the condensed-matter physics, the coupled cluster method has up to now been successfully applied to ultracold 
gases of bosonic atoms in traps [73, 74] and allowed to describe correlations beyond the mean-field regime in a Bose- 
Einstein condensate of thousands of atoms. However, the ground state of bosonic systems and bosonic many-body 
wave functions have much simpler form than the fermionic counterparts, and therefore the advantage of the CC 
method was not fully pronounced in that case. The CC approach has also found numerous applications in studies of 
spin-1/2 lattice models, both in one-dimensional chains and in two-dimensional square lattices (see, e.g.. Refs. [75, 76]). 
In the following, we will show that the CC method proves to be ideally suited to study the problem of many fermionic 
atoms in one-dimensional traps. 

The plan of our paper is as follows. Section II describes the theoretical framework, including the definition of 
the many-body Hamiltonian in subsection HA, and the discussion of the exact solution for the two-body case, and 
convergence of the energy with the size of basis set in subsection HB. This is followed by the summary of employed 
many-body approaches in subsection HC. Section HI presents the results on the convergence with the size of the 
one-particle basis set in subsection HI A, the convergence with the excitation level included in the coupled cluster 
Ansatz for the wave function in subsection HIB, the density profiles in subsection HIC, and the comparison with 
experiments in subsection HID. Section IV summarizes our paper and discusses possible future applications of the 
developed approach and further developments. Details of the analytic solutions of the two-body case of subsection H B, 
and the derivations of the corresponding convergence laws with size of the basis set are presented in appendices A-C. 


II. THEORETICAL FRAMEWORK 


A. The many-body Hamiltonian 


The Hamiltonian describing a system of N spin-I/2 fermions (atoms) in a one-dimensional harmonic trap reads 
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where Xi represents the coordinate of the i-th atom, m is the atom mass, uj is the frequency of the trap, and g is the 
strength of the two-body contact interaction. Throughout the paper we use units of energy and interaction strength 
that correspond to uj = m = h = 1. This amounts to measuring energies E in units of huj, lengths in units of the 
harmonic oscillator characteristic length aho = \/?i/(mw), and the interaction strength g in units of fitcoho- Obviously, 
the Hamiltonian is symmetric with respect to the transposition of two arbitrary space coordinates, implying that the 
eigenfunctions must transform according to an irreducible representation of the permutation group Sjsj. In this paper 
we limit ourselves to the spin-I/2 atoms, so that the wave function will be fully antisymmetric with respect to the 
simultaneous transposition of any space and spin coordinates of two atoms. From now on we assume that N = A^-l-iyj, 
describes the total number of atoms {N'f atoms with the spin projection 1/2, and atom with the spin projection 
— 1/2). We will consider either systems with Wj- = iVj,, that is, in the simplest spin singlet 5 = 0 state, or ones with 
7 ^ iVj, such that the total spin S = \N^ — Aj,], that is, in a spin-stretched (high-spin) configuration. 

In almost all many-body theories it is customary to use the so-called algebraic approximation [77], i.e. to use a 
finite set of one-particle functions. These one-particle functions are usually expanded in terms of some known basis 
functions. In our case the most natural choice of the basis functions is obviously given by the eigenfunctions of the 
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one-dimensional harmonic oscillator: 
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where Hn{x) are the Hermite polynomials. The integrals of the one-particle part of the Hamiltonian are diagonal, and 
given by the eigenenergies of the harmonic oscillator, while the two-particle integrals may be calculated numerically 
using the exact Gauss-Hermite quadratures. 

Obviously, eigenfunctions of the harmonic oscillator provide the basis set that is the simplest for numerical applica¬ 
tions. There are many other choices, that require more numerical effort, but should assure better numerical precision 
and convergence: an obvious example in the case of non-harmonic potential is to use the corresponding orbitals that 
are exact eigenfunctions of the one-particle Hamiltonian. Explorations of these choices go, however, beyond the scope 
of this paper. 


B. The two-body case and the convergence with the size of the one-particle basis set 


To investigate in detail the properties of the considered system it is useful to solve analytically the two-body case 
first. The exact solution of the two-body problem in a one-dimensional harmonic trap was found by Busch et al. [78] 
and by Franke-Arnold et al. [79]. For convenience and to introduce the notation we decided to include a detailed 
description of the solution of the aforementioned Hamiltonian in appendix A. The ground state wave function is 

«'(X, x) = C„C, Ui-e/2, 1/2, x^), (3) 


where X = (xi -I- Xi)l'j2 and x = (xi — X2)l\/2 are the center of mass and relative coordinates, and C„ and 
are the normalization constants for the X- and x-dependent portions. The corresponding total energy of the system 
(including explicitly the zero-point energies of the two particles) is then E = e + n + 1, while U{a, b, x) stands for the 
Tricomi function [80]. The energy of the relative motion e is determined from 
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where r(.) is the Gamma function. See appendix A for more details. 

To address the problem of the energy convergence with increasing number U}, of one-particle functions in the basis 
set, we have pursued two complementary strategies. The first follows the one developed by Hill [81] for the case of the 
helium atom, where the exact wavefunction, expanded in terms of an infinite sum over basis functions, is truncated 
at a fixed number Uf, of basis set functions. This leads to an approximated energy Cn^ converging to the exact value 
e as (see appendix B for details) 
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Note that the formula found by Hill, although derived explicitly for the two-electron case, has been successfully applied 
across a wide range of atomic and molecular calculations (see e.g. Ref. [82, 83]). Our aim here is similar, in that 
the above functional form of the two-body extrapolation formula will prove crucial for obtaining accurate results for 
systems with N > 2, as discussed below. 

Since the approach outlined above does not allow for the optimization of the coefficients in a smaller Hilbert space, 
we also followed a second strategy, based on an actual variational minimization of the energy in the space spanned 
by the states contained in the basis set. The result reads (see appendix G for details) 
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where x = g~^{e), as given by Eq. (4), x' = 9e[g“^(e)], and A = —1/{\/2'itx'). 

Note that both approaches yield the same functional form of the convergence formula, differing in coefficients that 
are specific to the two-body case, and leading terms in Eq. (5) and Eq. (6) are numerically confirmed to be the same. 
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C. Summary of the many-body approaches 


As already stated above, almost all many-body theories are based on the algebraic approximation [77], i.e., on the 
parametrization of the wave function by expansion in a finite set of basis functions. Since we deal with states with 
total spin S' = 0 or high-spin (spin stretched) states with the total spin S = \N-f — A(|^|, the first approximation to the 
exact wave function will be the Slater determinant built of one-particle functions obtained by solving the Hartree- 
Fock equations of the mean-field theory (cf. [84]). We will denote the N functions occupied in the reference Slater 
determinant by while the remaining one-particle functions that are not occupied in the reference determinant 

will be denoted by {((>p}p" jv+ii where ni, is the number of the harmonic oscillator functions used to expand the one- 
particle solutions of the Idartree-Fock equations. Note that the total number of one-particle functions is 2nb instead of 
rifo since 4>a and (j)p include the dependence on the spin coordinate through the spin functions. Obviously, the relation 
rib ^ N must hold. There is a theorem stating that if a set of one-particle functions spans the one-particle Hilbert 
space L^(R^), then the set of all A^-particle determinantal wave functions constructed from this one-particle set will 
span the antisymmetric part of the N particle Hilbert space [85]. In other words, the N particle wave function can 
rigorously be expanded in terms of the Slater determinants built from the complete set of one-particle functions. The 
same theorem holds in the algebraic approximation, and it is the basis of exact diagonalization, otherwise referred to 
as the full configuration interaction (FCI) method. The difference between the mean-held and FCI energies is called 
the correlation energy. 

In the full conhguration interaction method the wave function of the many-fermion system is represented as 

^f=(I-HC)$, (7) 

where the Cl operator C may be written as a linear combination of /-particle excitation operators. 
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with coefficients Cp^ ;p^‘. Summation over the repeated lower and upper indices is assumed. Note that the indices a 
and p refer to the one-particle functions that are occupied or empty in the reference Slater determinant $, so Eq. (7) 
can be viewed as an operator representation of the expansion of the exact wave function in terms of the all possible 
Wparticle determinantal wave functions constructed from the {4>a}a=i®{4>pY'f!^N+i one-particle functions. 

The /-particle excitation operator is given by the product of one-particle excitation operators . The 

one-particle excitation operators are in turn dehned in terms of the conventional fermionic creation and annihilation 
operators e]) = aj,aa, where the creation and annihilation operators are defined with the respect to the physical 
vacuum |0), but applied to the Fermi vacuum 4> [86]. Note that the one-particle functions (pa and (pp are always 
orthogonal, so that the excitation operators commute. Note also that our approach is limited to the N particle 
Hilbert space (layer of the Fock space with the fixed number of spin-1/2 fermions). This means that the algebra of 
the excitation operators is not only commutative, but also nilpotent. The latter property easily follows from the fact 
that we can replace at most N one-particle functions pa hr the reference Slater determinant, and any further action 
will give zero as the result. An important consequence of this property is that any transcendental function of the 
excitation operators that has a well defined Taylor expansion will reduce to a finite polynomial. The Cl coefficients 
Cpp'"pP and the energy are obtained by diagonalizing the Hamiltonian matrix constructed from the matrix elements 
Him = . This explains why the FCI method is also referred to as the exact diagonalization 

method. 

The computational cost of the FCI method is prohibitive, as it scales as A^^(2nt,)^+^ for N ^ rib [87], which limits 
this method to small systems (small N) or small number of basis functions rib. One possible way to cure this problem 
is to limit the number of excitations included in Eq. (8), e.g., to single and double excitations. Such a truncation may 
lead to considerable savings of computer time, but unfortunately has a serious drawback. Any truncated Cl method 
is no longer size-extensive, which means the energy of two non-interacting systems is not the sum of the energies of 
these two systems. 

To overcome the size-extensivity problem of the limited Cl expansions, the coupled cluster (CC) method was 
introduced, first in the nuclear physics [58] and slightly later in quantum chemistry [59, 60] and in the electron gas 
theory [88]. In this method the wave function is given by the following exponential Ansatz 

4' = e^$, (9) 

where the cluster operator T is given by 
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Note that by comparison with Eq. (7) the following relation between the FCI and cluster operator must hold, f = 
ln(l + C), but since the algebra of the excitation operators is nilpotent, the logarithm, as well as the exponential 
function in Eq. (9), reduce to finite polynomials. The CC method is non-variational, in the sense that it does not 
involve any optimization over a set of variational parameters. The energy is given by the expression 

£: = , ( 11 ) 

while the cluster amplitudes are obtained by solving the following system of nonlinear algebraic equations 

0 = , / = 1,..., N. (12) 

Note that by virtue of the Baker-Campbell-Hausdorff formula the exponential factor reduces to multiple 

commutators of H and T, and one can prove that this commutator expansion is finite and contains at most four-fold 
commutators if only two-particle interactions are present in the Hamiltonian [58]. 

Obviously, the coupled cluster method including all excitations is fully equivalent to the FCI method, and its 
computational cost is as high. However, the truncated CC methods are much less time consuming, and unlike the 
truncated Cl methods, remain size-extensive. Due to the exponential form of the Ansatz, Eq. (9), the CC method 
truncated to single and double excitations effectively includes triply, quadruply, and higher excited determinants 
through the products, e.g. T 1 T 2 or . In the present work we will use the series of approximations with the cluster 
operator limited to single and double (CCSD) [89], single, double, and triple (CCSDT) [90], and single, double, 
triple, and quadruple excitations (CCSDTQ) [91-93]. The computational cost of these methods scales as n®, n®, 
and respectively. For a system consisting of a dozen of atoms methods including triple excitations cannot 
reasonably be used, while the experience gained for many-electron systems suggests that the triple excitations have 
an important contribution to the correlation energy. To reduce the computational cost of the CC method including 
triple excitations, the CCSD(T) method with the n\ scaling was introduced [94]. In this method the CCSD equations 
are solved iteratively, and the contribution of the triple excitations to the correlation energy is evaluated from the 
expression based on the many-body perturbation theory. It turned out that the CCSD(T) method is very accurate 
for many properties of atoms and molecules and, as for now, it is considered as the golden standard of quantum 
chemistry. As we show in the next section, this method can also be applied with success to the system of ultracold 
atoms in a one-dimensional harmonic trap interacting via a short-range contact type potential. 

In order to introduce the perturbation expansion of the correlation energy, often referred to as many-body perturba¬ 
tion theory (MBPT) [63], or Mpller-Plesset perturbation theory (MPPT) [95-97], it is useful to rewrite the equations 
for the cluster operator T in the following operator form [98, 99] 

f = h(e-*We^) = n{w+[W,f] + i[IT,f ]2 + i[tT,f ]3 + , (13) 

where W is the correlation (fluctuation) potential in the Mpller-Plesset partitioning of the Hamiltonian 

H = F + W, (14) 

F is the Fock operator, and [IF,T]„ denotes the n-fold nested commutator: [IT,T]o = W, \W,T]n+i = [[l^iT]. 
The nested commutator expansion in Eq. (13) terminates after the quadruple commutators since in our case the 
operator W contains only two-particle interactions. The resolvent superoperator 7Z is defined for arbitrary operator 
A as [98, 99] 
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and Ck denotes the one-particle energy associated with the one-particle function labeled by the index k. We assume 
that the energy of the highest occupied one-particle function is smaller than the energy of the lowest unoccupied in 
the reference determinant 4), so the superoperators IZn are always well defined. Note that for a given X the operator 
Y = TZniX) can be viewed as a formal solution of the equation 


{eP^\-::'LimF,y] + xm = o- 


( 17 ) 
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Obviously, this solution is unique if we assume that Y belongs to the linear span of the excitation operators ■ 

The many-body perturbation expansion of the energy is obtained by introducing the following parametrization of 
the correlation operator IT, 

W XW, (18) 

where the complex parameter A is introduced to derive the perturbation expansions of Eqs. (11) and (13). Obviously, 
the physical value of A is equal to one. By substituting the parametrized Eq. (18) into Eq. (14) the cluster operator 
became dependent on A, and can be expanded as a power series in A 

OO 

f = (19) 

fe=i 

Substituting the above expansion in the expression for the energy (11) and collecting all terms at A" gives the 
expression for the nth-order correction to the energy in the many-body perturbation theory. The cluster operators 
necessary to evaluate this expression are obtained by substituting Eq. (19) and again collecting all terms at A". 

It follows immediately from this short sketch of the derivation that the MBPT and CC theories are strongly 
connected. The analysis reported in Ref. [99] shows that the CCSD, CCSDT, and CCSDTQ are valid through the 
third, fourth, and fifth order of the many-body perturbation theory. The golden standard of quantum chemistry, the 
CCSD(T) method, is also valid through the fourth order of MBPT, although it is computationally less demanding 
than the full CCSDT theory. 

One of the observables that can be obtained with the presented method is the density profile of the trapped atoms. 
At a point xq it may be obtained as the expectation value of the operator ~ *o) within the FCI calculations, 

and from the Hellmann-Feynman theorem as the first derivative of the energy in the presence of the perturbation 
given by the same operator in the CC calculations. In practice, we compute this derivative by using the finite-field 
approach with a perturbation of the form (j)i[xo)4>jixo) added to the *j-th element of the one-particle Hamiltonian 
matrix. 

The full configuration interaction and coupled cluster calculations were performed with the customized versions of 
the HECTOR [100] and ACESH codes [lOI], respectively. The CC and truncated Cl calculations use the reference 
state $ built of orbitals obtained with the restricted Hartree-Fock method for the case, whereas orbitals 

obtained with the unrestricted Hartree-Fock method are employed for the ^ case. 


III. NUMERICAL RESULTS AND DISCUSSION 
A. Convergence with the size of the one-particle basis set 

Examples of a slow convergence of the results with respect to the basis set size are well-known in the literature. 
The conventional exact diagonalization calculations for solving the electronic Schrddinger equation converge as L~^, 
where L is the highest angular momentum present in the one-electron basis set [81, 102]. Even slower convergence 
{L~^) was observed for some relativistic corrections arising from the perturbative approach based on the Breit-Pauli 
Hamiltonian [103]. 

In order to check that the convergence formulas of Eq. (5) and Eq. (6) for the two-body problem hold, we compare 
their predictions with the exact results in Fig. 1. Two sets of calculated energies relative to the exact value are 
presented as a function of the one-particle basis set size: non-variational energies obtained as the expectation value with 
the truncated exact wavefunction [Eq. (B5) in appendix B] and variational energies equivalent to exact diagonalization 
results. The former one should be described by the asymptotic expansion of Eq. (5), whereas the latter ones should 
coincide with the asymptotic expansion of Eq. (6). The asymptotic formulas with increasing number of terms are 
presented. The first terms of the asymptotic expansions give a reasonable estimate of the convergence rate and upon 
inclusion of the second and third terms the differences between the analytical formulas and calculated values are 
greatly reduced. Therefore, the formulas (5) and (6) are valid and can serve as a guide for further investigations. 

To evaluate the adequateness of the functional formula given by Eqs. (5) and (6) for extrapolation of the variational 
energies to the complete basis set limit, Fig. 1 presents also the fit of the variational results to the formula Eoo + 
Arifj -I- This extrapolation formula behaves extremely well [see also the inset of Fig. 1]. For 

instance, the extrapolated energies in the complete basis set limit Eoo for g = 1 and g = h are 1.306744 and 1.726780, 
whilst the exact values equal to 1.306745 and 1.726771, respectively. Note that in the biggest basis set available, 
rifo = 200, the calculated energies are 1.310545 and 1.745325, respectively, so that the accuracy gain of three order of 
magnitude is impressive and very important. 
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FIG. 1. Absolute error with respect to the exact result of the ground state energy of the 1+1 system for g = 5. Comparison of 
the non-variational energy (blue squares) obtained as the expectation value with the truncated exact wavefunction [Eq. (B5) 
in appendix B] and variational energy (red dots) equivalent to exact diagonalization results with the asymptotic formulas of 
Eqs. (5) and (6). Results with inclusion of the leading-order term only are given as the short-dashed lines, results with inclusion 
of the first two terms are given as the dashed lines, and with inclusion of all three terms are given as the dot-dashed lines. 
The solid black line is the fit of the variational results to the formula, E^o + and the inset shows the 

performance of this ht over full set of data. 



9=5 


9=20 


0.00 0.05 0.10 0.15 

l/n^^ 


0.20 



FIG. 2. Absolute error with respect to the complete basis set limit of the ground state energy of the 2+2 (a) and 3+3 (b) 
systems as a function of the one-particle basis set size for several interaction strengths obtained with the FCI method. Lines 
correspond to fits to the extrapolation function given by Eq. (20). 


Given the very slow convergence rate of the two-body energy with the number of the basis functions ri},, it is very 
important to consider this convergence in the many-body case. One can expect that the convergence pattern obtained 
for the two-body case is almost certainly valid for the many-body case. This fact was observed in the electronic 
structure calculations on atoms and molecules, and is related to the fact that the correlated pair functions of two 
electrons have the same behavior around the electrons coalescence points as in the discussed two-body example. 
Clearly, this analytic behavior determines the convergence rate of the results towards the exact value. A very similar 
situation is found for the partial wave expansion of the exact solution of the Schrodinger equation. The well-known 
L~^ convergence pattern has rigorously been derived only for the helium atom [81, 102], but was successfully applied 
also for many-electron atoms/molecules [104, 105] and is widely accepted to be universal (c.f. the discussion given by 
King [106]). 

Following the discussion above and guided by the expressions (5) and (6), we decided to adopt the following 
















FIG. 3. Energy spectra (corrected by the ground-state energy of the noninteracting system Ep) for the 2+2 (a) and 3+1 (b) 
systems as a function of the interaction strength g obtained with the FCI method in the basis of rit = 20 one-particle functions 
and in the complete basis set limit nj, = oo. 


three-term extrapolation formula 


E-ni, ~ Ec 
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where E^c is the extrapolated energy to the infinite number of the basis functions, En^ is the energy computed with 
rib basis functions, and A and B are the fit parameters. To check the correctness of the above expression we performed 
FCI calculations for the ground state of a balanced two-component Fermi gas. We have found that the convergence 
pattern with the number of the basis functions exactly follows our extrapolation formula. The accuracy of our formula 
is illustrated in Fig. 2. There we show the absolute error with respect to the complete basis set limit in the ground 
state energy for the 2+2 and 3+3 systems (A^-i- + means atoms with the spin projection 1/2, and Ni atom 
with the spin projection —1/2) as a function of the number of one-particle basis functions for several interaction 
strengths obtained with the FCI method, as well as the error plotted as a function of rib according to Eq. (20). The 
agreement between the computed points and the analytical fits to the extrapolation function is excellent, supporting 
the correctness of our extrapolation scheme. 

In order to get saturated results, careful studies of the convergence are necessary. As stated above, the system of 
many identical fermions with spin-1/2 in a ID trap interacting via the contact potential has an even worse convergence 
than the one observed for the electronic Schrodinger equation [81, 102, 103], and apparent convergence may be 
observed. This is illustrated in Fig. 3 where we show the energy spectra for a few atoms in a trap obtained with a 
small number of the basis functions {rib = 20 as in Ref. [48]) and in the complete basis set limit. An inspection of 
this figure shows that the results obtained with 20 basis functions are far from the complete basis set limit obtained 
from Eq. (20), especially in the limit of intermediate and large interaction strength. The apparent convergence of 

_ 1 /Q 

the results observed by the Authors of Ref. [48] is solely due to the pathological convergence pattern as . As 
will be shown in the next section, a quantitative picture of the physics, and a quantitative agreement with various 
experimental data, can only be obtained if the extrapolation to the complete basis set limit is properly done. 

As seen in Eigs. 2 and 3, both the absolute and relative errors due to the finite basis set size are increasing with the 
increasing interaction strength g. In the limit of very strong interactions, particularly approaching the unitary limit 
of Tonks-Girardeau (TG) gas (when the interaction strength is infinite, g = oo), this error can become larger than the 
separation between the ground and excited states and will lead to artefacts in the description of physical phenomena 
both at zero and hnite temperatures. The Authors of Ref. [48] investigated with the exact diagonalization method 
the few-fermion physics in the TG limit when the spectrum is quasi-degenerate, but the careful convergence analysis 
reveals that the results obtained with the basis set as small as rib = 20 give only a qualitative picture in this regime 
of interaction strength. 

Figure 4 shows the energy spectrum for the example of the 3+1 system as a function oil/ g obtained with the FCI 
method in the limit of the strong interaction in the basis of rib = 20 one-particle functions and in the complete basis set 
limit rib = oo. The complete basis set limit results were obtained by extrapolation from the numerical results in basis 
of 30, 40, 50, and 60 functions. In the TG limit the four states become degenerate. Unfortunately, any calculation 
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FIG. 4. Energy spectrum for the 3+1 system in the limit of the strong interaction as a function of l/g obtained with the FCI 
method in the basis of nt, = 20 one-particle functions and in the complete basis set limit nt = oo. The dashed black lines are 
the linear fits to the points in the complete basis set limit and the solid red lines are the analytilcal results obtained with the 
perturbation theory (PT) in Ref. [53]. 



FIG. 5. Gonvergence of the critical coupling l/g^ (at which states in the ground manifold become degenerate) as a function 
of the one-particle basis set size obtained with the FGI method. The various symbols indicate systems with N = 
particles and the solid line is given by the first term of Eq. (G5) in appendix G. 


with a finite number of one-particle basis set functions will always locate a crossing of these states at a finite g, 
potentially leading to incorrect conclusions. For rib = 20 the crossing appears at g = 17.1 and corresponds to the 
large error of 0.4 a; in the ground-state energy. Therefore, the extrapolation to the complete basis set limit is not only 
needed to improve accuracy, but is necessary to have a quantitatively correct picture of the physics. The numerical 
results in the complete basis set limit agree perfectly well with the analytical results obtained in the vicinity of the 
Tonks-Girardeau limit with the perturbation theory [50, 53]. The slopes of the curves fitted to the numerical data 
are {—7.14,—3.58,—1.18,0}, in a very good agreement with the exact contact coefficients {—7.08,-3.58,-1.18,0}. 
This confirms that the presented extrapolation scheme works also very well in the regime of strong interactions 
allowing to approach the TG limit with the accurate finite basis set calculations. Interestingly, the contact coefficients 
corresponding to non-converged results, even in basis set as small as Ub = 20, are surprisingly close to the correct 
values for the TG limit. 

As we have shown above, the energy spectrum in the TG limit obtained in the calculations with a single one-particle 








10 



g/hcoa^^ 


FIG. 6. Correlation energy for the 2+2 (a) and 3+3 (b) systems as a function of the interaction strength obtained with the 
FCI method for several one-particle basis set sizes and complete basis set limit. 


basis set has always an artificial crossing of the states at a finite value of the interaction strength g. In Fig. 5 we 
plot the values oi 1/g for which this crossing occurs as a function of the one-particle basis set size for a few systems. 
Interestingly, the interaction strength at which the states in the ground manifold become degenerate depends solely 
on the number of the basis functions and not on the number of atoms, neither on their state. This observation agrees 
with our prediction on the convergence of 1/g for a given energy (in this case the energy of the ferromagnetic state) 
with the size of the one-particle basis set to be independent in the leading order on both interaction strength and 
energy. The leading term of the convergence formula for the two-body problem, derived analytically in appendix C, 
is shown as a solid black line in Fig. 5. This observation suggests a second approach to get complete basis set limit 
results and to reach the TG limit with accurate finite basis set calculations, that is, instead of extrapolating the 
energy, one can extrapolate 1/g for fixed energies by using the universal convergence formula. Exactly the same 
convergence coefficient valid for all presented few-body cases and observed universality suggest even a third approach 
to get accurately converged results with finite basis set calculations at a given interaction strength g. The idea is to 
use in the numerical calculations a renormalized coupling constant gm explicitly dependent on the high-energy cut-off 
set by rib, which reproduces exact two-body result in a smaller Hilbert space spanned by nt basis functions. In this 
way, the inaccurate description of the short-distance two-body physics introduced by the high-energy cutoff may be 
cured order by order. A similar idea was discussed, e.g., in Refs. [107, 108]. 


B. Convergence with the excitation level included in the coupled cluster Ansatz for the wave function 

Having solved the problem with the convergence of the results with the number of basis functions, we now turn to 
the effect of the truncation of the excitation in the cluster operator, Eq. (10), on the correlation and total energies. 
We start the discussion of our results with Fig. 6 where we report the correlation energy for the 2+2 and 3+3 systems 
as a function of the interaction strength for two one-particle basis set sizes and the complete basis set limit. An 
inspection of Fig. 6 shows that the basis set dependence of the energy is relatively weak in the weakly correlated 
repulsive regime {g > 0), and very pronounced in the strongly correlated attractive case (g < 0). These significantly 
different behaviors result from the fact that in the limit of the strong repulsion even distinguishable fermions tend to 
occupy different one-particle states and the total wave function approaches the structure of the ferromagnetic state 
(the Tonks-Girardeau gas limit). On the other hand, the fermions with the opposite spin projection tend to pair in 
the case of the attractive interaction and become tightly bound (hard-core) bosonic dimers in the limit of the strong 
attraction (the Lieb-Liniger gas limit). The description of these tightly bound pairs is challenging for calculations in 
the one-particle functions basis sets and explicitly correlated methods could potentially overcome the problem and 
signihcantly accelerate the convergence. It is important to stress at this point that the mean-field energy converges 
very fast with the number of one-particle functions, so the convergence problems are related to the basis saturation 
of the correlation energy. 

The correlation energy presented in Fig. 6 diverges linearly with the interaction strength g for both negative 
and positive values. One should note that in the weak interaction regime the distinction between the mean-field and 
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FIG. 7. The percentage of the ground state correlation energy reproduced with the many-body perturbation theory MBPT2, 
MBPT4, and the truncated configuration interaction method CISD in the 2+1 (a), 2+2 (b), 3+2 (c), and 3+3 (d) systems as 
a function of the interaction strength. Values are obtained for a selection of basis set sizes and extrapolated to the complete 
basis set limit. 
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FIG. 8. The percentage of the ground state correlation energy reproduced at the CGSD, CCSD(T), CGSDT, and CGSDTQ 
levels of the coupled cluster theory in the 2+1 (a), 2+2 (b), 3+2 (c), and 3+3 (d) systems as a function of the interaction 
strength. Values are obtained for a selection of basis set sizes and extrapolated to the complete basis set limit. 


correlation energies is reasonable. For intermediate and especially strong interaction regimes the mean-field description 
fails completely and diverging correlation energy compensates unphysically diverging mean-field energy, so no physical 
meaning should be attributed to this divergence. However, this behavior of the mean-field and correlation energies 
affects the performance of post-mean-field methods in the strong interaction regime (including the standard coupled 
cluster method), which start from the mean-field solution and tend to recover correlation energy. 

We now turn to the applicability of the many-body perturbation theory and truncated configuration interaction 
expansions for the energy calculations. In Fig. 7 we report the percentage of the ground state correlation energy 
reproduced with the second-order and fourth-order many-body perturbation theory, MBPT2 and MBPT4, and the 
configuration interaction method limited to single and double excitations, CISD, for the 2+1, 2+2, 3+2, and 3+3 
systems as functions of the interaction strength. As expected from the electronic structure calculations on atoms and 
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TABLE I. The percentage of the ground state correlation energy reproduced at the CCSD, CCSD(T), CCSDT, and CCSDTQ 
levels of the coupled cluster theory and the many-body perturbation theory MBPT2, MBPT4, and the truncated configuration 
interaction methods CISD in the 2+1, 2+2, 3+2, and 3+3 systems for the different interaction strengths. N/D (no data) 
indicates that calculation was not feasible with used software. 


system 

CCSD 

CCSD(T) 

CCSDT 

CCSDTQ 

MBPT2 

MBPT4 

CISD 




g = 

2 




A = 2 + 1 

99.60 

100.07 

100 

- 

117.33 

100.31 

99.42 

A = 2 + 2 

99.30 

100.10 

99.95 

100 

115.26 

101.49 

98.70 

A = 3 + 2 

99.35 

100.07 

99.97 

N/D 

113.48 

100.18 

97.40 

A = 3 + 3 

99.23 

100.06 

99.94 

99.999 

120.07 

101.65 

96.17 




g = 

4 




A = 2 + 1 

97.82 

99.69 

100 

- 

92.42 

98.28 

95.31 

A = 2 + 2 

98.53 

100.72 

99.90 

100 

153.09 

127.21 

95.39 

A = 3 + 2 

96.89 

100.18 

99.86 

N/D 

97.54 

95.05 

87.18 

A = 3 + 3 

98.14 

100.56 

99.84 

99.997 

144.33 

116.83 

90.64 




g = 

-4 




A = 2 + 1 

98.33 

99.46 

100 

- 

76.51 

100.08 

98.24 

A = 2 + 2 

99.53 

101.97 

102.64 

100 

68.49 

96.08 

87.60 

A = 3 + 2 

97.96 

100.10 

100.35 

N/D 

97.54 

95.05 

87.18 

A = 3 + 3 

99.79 

102.35 

102.99 

100.003 

144.33 

116.83 

90.64 


molecules the performance of the MBPT2 and MBPT4 methods is not very good, and quite erratic as a function 
of the interaction strength. Given the fact that for atoms and molecules, the many-body perturbation theory is 
divergent [109-111] as perturbation theories applied in a different physical context [112, 113], it is not surprising that 
the MBPT4 theory is not a big improvement over the MBPT2 approach, despite a much higher theoretical complexity 
and computational time scaling with the size of the basis, for MBPT2 vs. nl for MBPT4. Finally, we note that 
the variational CISD results do not offer a good advantage over the MBPT results. 

Since the MBPT and limited Cl theory perform badly, one may ask if a selective infinite-order summation of some 
MBPT diagrams with different variants of the coupled cluster theory will improve the situation. This is indeed the 
case. Fig. 8 shows the percentage of the ground state correlation energy reproduced with the CCSD, CCSD(T), 
CCSDT, and CCSDTQ methods for the 2+1, 2+2, 3+2, and 3+3 systems as functions of the interaction strength. 
First of all note that for the 2+1 system CCSDT is equivalent to the FCI theory, and for the 2+2 system the same 
statement is valid for CCSDTQ. This means that for these systems these particular methods will reproduce 100% 
of the energy. It may be surprising at first glance that some of the approximate variants of the coupled cluster 
theory overestimate the total energy of the system. This is due to the fact that the truncated CC methods are not 
variational, and one can obtain e.g. 101% of the energy. An inspection of Fig. 8 shows that all CC methods perform 
very well, although there is no clear convergence pattern with the excitation operators included in the CC wave 
function. Indeed, the CCSD method tends to slightly underestimate the energy, while methods including triple (and 
possibly quadrupole) excitations tend to slightly overestimate. In general, the percentage errors are of the order of a 
few percent. Note also that the CCSDT method does not offer big advantage over the CCSD(T) approach, despite of 
much higher computational requirements. This is in line with the results obtained for atoms and molecules, see for 
instance Refs. [103, 114]. Comparison of Figs. 7 and 8 shows that the cluster expansion of the wave function and the 
summation of the diagrams involving the single and double excitation in the MBPT theory to infinite order is crucial 
for an accurate description of many-atom systems in a one-dimensional harmonic trap. 

The discussion of the last two paragraphs is well summarized in Table I, where we show the CC, MBPT, and CISD 
results for selected values of the interaction strength. An analysis of the numerical results reported in this table 
confirms a very good performance of various variants of the coupled cluster theory, as opposed to the erratic behavior 
of the many-body perturbation theory, and the poor performance of the configuration interaction method limited to 
single and double excitations. 

Finally, we note that the performance discussed above of the truncated CC methods shows a weak dependence 
on the number of one-particle basis functions. Similarly, we do not observe any increase of the relative error with 
the number of atoms, although our comparison with the FCI results is restricted to the maximum of six atoms. For 
a larger number of atoms we can evaluate the contributions of the noniterative triple excitations in the CCSD(T) 
method. Figure 9(a) shows the interaction energy in the ground state of the 10+10, 20+20, 30+30, and 40+40 
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FIG. 9. (a): The ground state energy (relative to the ground-state energy of the noninteracting system Ep) of the 10-1-10, 
20-1-20, 30-1-30, and 40-1-40 systems as a function of the interaction strength obtained at the CCSD and CCSD(T) levels of the 
coupled cluster theory extrapolated to the complete basis set limit, (b): The percentage of the ground state interaction energy 
accounted by the connected triples contribution calculated non-iteratively within many-body perturbation theory as a function 
of the interaction strength for the same systems. 


systems obtained at the CCSD and CCSD(T) levels of the coupled cluster theory. Figure 9(b) presents the percentage 
of the ground state interaction energy accounted by the difference between energy obtained with the CCSD(T) and 
CCSD methods, i.e. by the connected triples contribution calculated noniteratively within the many-body perturbation 
theory. Interestingly, for a given interaction strength the contribution of the excitations higher than double to the 
ground-state interaction energy is becoming less important with increasing number of atoms in the system. Based on 
this fact we predict that the performance of the CCSD(T) method for larger systems is as good as for the investigated 
small systems and the lack of higher excitations in the coupled cluster wave function does not introduce any significant 
errors. 

The energy spectrum of the investigated systems for strong repulsive interactions becomes quasi-degenerate 
(cf. Figs. 3 and 4) and the correlation energy diverges for both strongly attractive and repulsive interactions (cf. Fig. 6). 
These two effects lead to the convergence problem of the standard single-reference coupled cluster method starting 
from the antiferromagnetic reference state as used in the present study, restricting the interaction strength g that can 
effectively be used in actual calculations to intermediate values between -6 and 6. The possible solution to overcome 
this problem is the use of the multireference version of the coupled cluster theory or the ferromagnetic reference state 
for calculations in the limit of the strong repulsive interaction. 


C. Density profiles 

The density profiles of the trapped clouds are important observables that can be measured experimentally and used 
to monitor the evolution of the interatomic interactions and resulting states of a many-body system. 

In Ref. [II5] we have analyzed the density profiles for the 3-1-3 system with the FCI method and for the 15-1-15 
system with the CC approach. We have shown that both methods perfectly reproduce the density profiles of the 
non-interacting gas, and that the FCI calculations allow to approach both the limit of strong attraction (the Lieb- 
Liniger gas) when the fermionic atoms of opposite spin projection pair into hard-core bosonic dimers, and the limit 
of strong repulsion (the Tonks-Girardeau gas) when even distinguishable fermions must occupy different one-particle 
levels. The CC method allows to describe the evolution of the density prohles in the range of intermediate values 
of the interaction strength. The densities obtained with two methods agree with predictions of the local density 
approximation applied to the solution of the Gaudin-Yang integral equations describing a homogeneous gas [38] 
providing a further confirmation of the validity of the local density approximation for investigating a trapped ID gas. 
Here, we present other examples for a few fermion systems obtained with the FCI method and for the many fermion 
systems calculated with the CC approach. 

Figure 10 shows the density profiles for the 2-|-2 and 3-1-1 systems obtained with the FCI method for three interactions 
strengths: strongly repulsive {g = 50), moderately repulsive {g = 5) and strongly attractive {g = —10), together with 
the analytic results. In the balanced case, it is straightforward to find analytical expressions for the limiting cases 
of strong attraction (nLL(a:) = the Lieb-Liniger gas of hard-core bosonic dimers at g = —oo. 
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FIG. 10. Density profiles of a two-component Fermi gas obtained with the FCI method for the 2+2 (a) and 3+1 (b) systems 
for the three interaction strengths g — 50, g = 5 and g = —10. The analytic results for the limiting cases of strong attraction 
{g = —oo), strong repulsion {g = oo), and no interaction {g = 0) are also presented. 




FIG. 11. Density profiles of a two-component Fermi gas obtained with the GCSD(T) method for the 5+5, 10+10, and 20+20 
systems and the interaction strength g = 1 (a), and for the 5+5 system with different interaction strengths (b). The analytic 
results for the limiting cases of strong attraction {g = —oo), strong repulsion {g — oo), and no interaction {g = 0) are also 
presented in panel (b). 


with (pi{x) the i-th eigenfunction of the harmonic oscillator for a dimer of mass 2m), strong repulsion ( 71 x 0 ( 2 ;) = 
Si=o^ |;pi(a;)p for the Tonks-Girardeau gas of ” ferinionized” fermions at g = 00 ), and no interaction (no(x) = 
\g}i{x)\'^ for g = 0). In the case of strongly repulsive interaction, the FCI results are indistinguishable from 
the analytic result confirming perfect performance of the FCI method in this limit of the interaction strength. The 
regime of the strong attraction is much harder to describe due to the presence of strong two-body correlations when the 
hard-core bosonic dimers are formed. The presented results for g = —10 approach the analytic results but obtaining 
fully converged results for much more attractive interaction strengths, even with the extrapolation used within the 
paper, is very challenging. A possible solution to overcome the very slow convergence in the limit of strong attractive 
interaction can be to use the explicitly correlated method or to include the formation of the hard-core bosonic dimers 
directly in the structure of the wave function. 

Figure 11(a) shows the density profiles for the 5+5, 10+10, and 20+20 systems obtained with the CCSD(T) method 
for the interaction strength g = 1. With an increasing number of atoms the evolution of the size of the cloud can 
be observed. As we have shown in Ref. [115], the overall shape of the profiles can be described by the typical 
Thomas-Fermi profile of an inverted parabola. The relative size of the density oscillations are smaller and smaller 
with the increasing number of atoms, and the density profile approaches the exact shape given by the local density 
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FIG. 12. (a): The ground state energy (relative to the ground-state energy of the noninteracting system Ef) of the N atoms 
consisting of a single impurity with an opposite spin interacting with an increasing number of identical fermions for various 
interaction strengths g obtained with the FCI method is compared to the measurements of Ref. [32]. (b): The difference 
between theoretical and experimental values with conservatively estimated error bars of numerical results, (c): The separation 
energy obtained with the FCI method and measured in experiment [33]. 


approximation in the thermodynamic limit. Figure 11(b) shows the density profiles for the 5-1-5 system obtained with 
the CCSD(T) method for the three interactions strengths g = —4, g = 2, and g = 6- The present CC calculations 
are restricted to intermediate values of the interactions strength, which however extend far beyond the mean field 
regime [115]. Moreover, the accessible range is sufficient to observe significant modifications of the shape and the size 
of the cloud. The evolution of the cloud’s shape towards the analytic results of the limiting cases of strong attraction 
{g = —oo) and strong repulsion {g = oo) is clearly visible in Fig. 11(b). Finally, the CC method allows one to address 
much larger systems than the FCI approach. 


D. Comparison with the available experimental data 

All the results reported thus far strongly suggest that the approximate variants of the coupled cluster method 
combined with the efficient extrapolation schemes based on the rigorous analytical formulas derived for the two-body 
case perform very well for a few-body and many-body systems of identical spin-1/2 atoms interacting via the short- 
range contact potential in a ID harmonic trap. However, the most stringent test of the accuracy of any theoretical 
model is the comparison with precision experiments. 

We report such a comparison in Fig. 12. We decided to compare the best of our results, i.e. the FCI results, but 
any variant of the coupled cluster method that would include triple excitation would lead to the same results within 
1 %, indistinguishable from the FCI values on the scale of the plot. Panel (a) of this figure shows the comparison 
between the measured [32] and computed energies (with respect to the ground-state energy of the noninteracting 
system Ep) of a system of = A^-|- -I- 1 atoms consisting of a single impurity interacting with a Fermi sea of N-f 
identical fermions for various repulsive interactions. An inspection of this figure shows that the agreement between 
theory and experiment is indeed excellent. Note parenthetically that we include on this hgure both the experimental 
error bars and the computational uncertainties A^; conservatively estimated from the extrapolated energies and the 
energies computed with the largest ni, = 200 basis functions, Ap = En ^=200 — A'oo- On the scale of this graph the 























16 


experimental and theoretical error bars are indistinguishable, so in the panel (b) of this figure we report the differences 
between the theoretical and experimental energies with the respective error bars. The agreement viewed in this way is 
also very good, and the energy difference always lies within the combined theoretical and experimental uncertainties. 

Less satisfactory is the agreement between theory and experiment for the separation energies AE{N) = fi{N) — 
/i*(7V), where n{N) = E{N) — E{N — 1) is the chemical potential, E{N) is the extrapolated energy for the N- 
body system in the weakly attractive regime, and /r*(iV) is the chemical potential of the noniteracting system. This 
is illustrated in the panel (c) of Fig. 12, where we compare our theoretical results with the experimental data of 
Ref. [33]. An inspection of this figure shows that a relatively good agreement (but not within the experimental error 
bars) is observed for odd values of N, and important disagreement by a factor of roughly 3/2 for even N. Note that 
the experimental technique used in Ref. [33] is different than in Ref. [32]. The latter one is based on an accurate 
non-destructive measurement of the RF spectrum of the systems, whereas in the former one the system is probed by 
deforming the trapping potential and by observing the tunneling of particles out of the trap. We believe that the 
main source of disagreement in the second case comes from the perturbed character of the ID harmonic shape of 
the trap during measurements and an approximate determination of the interaction strength in the experiment [33]. 
Indeed, the results of the calculations with a smaller (in the absolute value) interaction strength g are in a much 
better agreement for both even and odd values of N, as it has also recently been shown by the authors of Ref. [54]. 


IV. SUMMARY AND CONCLUSIONS 

In this paper we have reported the first application of the ab initio methods of quantum chemistry to systems of 
many interacting fermions in a one-dimensional harmonic trap. Our results can be summarized as follows: 

• The behavior of the two-body energy as a function of the number of single-particle functions nj, included in 
the calculations for a fixed interaction strength has thoroughly been analyzed and an extrapolation formula has 
been derived. The convergence with rib is pathologically slow and is reported in the leading order as rij, ' . 

• The convergence of the interaction strength as a function of the number of single-particle functions nt for 
calculations at fixed energy in the two-body case have also been analyzed and an extrapolation formula which 
does not depend on the energy and on the number of particles in the lowest order has been derived. 

• The convergence with the number of basis set functions Ub in the many-body case has been analyzed and a 
suitable three-point extrapolation formula has been found. 

• The importance of the use of converged results to correctly describe the physics of the trapped Fermi systems 
has been pointed out. Shortcomings of the approaches lacking a proper analysis of the convergence issues have 
been shown. 

• Comparison of several quantum many-body methods has been reported. In particular, a careful consideration 
has been given to the level of the excitation present in the wavefunction in the coupled cluster method. An 
optimal method has been chosen and used throughout the rest of the study. 

• The limitations of the adopted coupled cluster implementation with the standard single-reference state have 
been discussed together with the analysis of the behavior of the decomposition of the total energy into the 
mean-field and correlation parts in the limits of strong repulsive and attractive interactions. 

• The coupled cluster method restricted to single, double, and noniterative triple excitations CCSD(T) has been 
applied to describe systems of up to 80 fermionic atoms. 

• Density profiles have been computed and analyzed in the full spectrum of the interaction parameter and the 
number of atoms. The transition between the Lieb-Liniger gas of hard-core bosonic dimers and the Tonks- 
Girardeau gas has been observed with the FCI method. 

• Comparison with the available experimental data including estimates of the computational uncertainties has 
been provided. Very good agreement between the theory and experiment has been pointed out for the precision 
measurements based on the RF spectroscopy whereas the perturbed character of the ID harmonic shape of 
the trap during measurements based on the observation of the tunneling of particles out of the trap has been 
confirmed by confrontation with our accurate ab initio calculations with estimated error bars. 

The presented numerical approach allows us to get an important insight into the ongoing experiments on the trapped 
fermions. By extending the number of atoms from less than ten within the exact diagonalization approach (see, e.g., 
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Ref. [56]) to many tens within the coupled cluster method, it is becoming possible to investigate the transition between 
few- and many-body systems. This is particularly important for a good understanding of the emergence of bulk matter 
properties, crucial across all areas of physics. 

Recently, we have applied the numerical approach developed here to investigate the properties of a two-component 
Fermi gas trapped in one dimension [115]. We have addressed the question of how the observables, such as the energy, 
the chemical potential, the pairing gap, and the density profile, evolve as the number of particles increases from very 
few to many tens. We have found that the energy converges surprisingly rapidly to the many-body result for all 
interaction strengths between minus and plus infinity, covering the whole range from the molecular bosonic Tonks 
gas to the atomic (fermionic) one. On the other hand, we observed the emergence of a non-analytic behavior of the 
pairing gap only when a substantially larger number of particles is present in the trap. 

We believe that the results presented here and in Ref. [115] on several fermions in a harmonic trap will be followed 
by many applications of the proposed numerical approach to investigate interesting physics in different systems, ge¬ 
ometries, and dimensions. We foresee that the coupled cluster method will describe atoms trapped in other potentials, 
such as double wells, arrays of microtraps, or even optical lattices equally well. The method should work also for 
more complex two-body interaction potentials, e.g., including van der Waals or long range 6dipole-dipole interactions. 
What is more, the method can be generalized to handle atoms trapped in two or three dimensions. The pathologically 
slow convergence with the number of single-particle basis functions included in the calculations can be overcome by 
using explicitly correlated basis functions in analogy with methods developed in quantum chemistry [116]. The use of 
the explicitly correlated methods can be crucial for an accurate description of systems in a higher number of dimen¬ 
sions. The range of interaction strengths which can be used in calculations with the single-reference coupled cluster 
method is restricted to intermediate values. This can be overcome by using a ferromagnetic (high-spin) reference 
state instead of the antiferromagnetic (low-spin) one which is typically used in the standard coupled cluster method. 
This would allow to describe many tens of interacting atoms in the strongly repulsive regime, and therefore would be 
ideal for describing atomic clouds in the Tonks-Girardeau limit. Finally, one can employ the coupled cluster method 
based on the high-spin reference state to describe many strongly-interacting bosonic atoms and investigate in detail 
the fermionization process in such a system. 
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Appendix A: The exact solution of the two-body case 


Let us consider the Hamiltonian (1) with = 1. By introducing new coordinates x = "^(a^i ~ ^ 2 ) (relative 

motion) and X = + X 2 ) (centre-of-mass motion) one arrives at the formula 





vI/(x,A)=0, 


(Al) 


where g = glV^- In these variables the exact wavefunction, 4'(a;, A), separates into a product of functions </> and ^ 
dependent only on x and X, respectively. These functions obey the following differential equations 

- + ]^x^4>{,x) + g5{x)4>{x) = 

/ \ (■^^) 

where primes denote the usual differentiation over the corresponding variables. The equation for the centre-of-mass 
motion can immediately be solved, as it coincides with the Schrodinger equation for the quantum harmonic oscillator. 

Further in the text we shall use the shorthand notation, -I- ^x'^ + gS{x) — 1/2, so that Hj;(j){x) = e<j){x). 

The Hamiltonian, is invariant with respect to the spatial inversion, i.e., x —> —x. Therefore, its eigenfunctions 
possess a definite parity (even or odd). Eigenfunctions that are of odd parity must vanish at a; = 0. 
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A substitution (t>{x) 


= 


f{x) in the first equation of (A2) gives the corresponding differential equation for fix) 


f"{x) - 2xf'{x) + 2e/(x) = 2gSix)fix), 


(A3) 


First, we solve the homogeneous differential equation, i.e., without the Dirac delta term on the right-hand-side, which 
is well-known and the general solution can be written as a linear combination of two functions 

Cu Ui-e/2, 1/2, x^) + Cm M{-e/2, 1/2, x^), (A4) 

where Cjj and Cm are some constants necessary to satisfy the initial conditions, and M and U are the Kummer and 
Tricomi hypergeometric functions [80] , respectively. Returning to the original (inhomogeneous) equation, we consider 
first the odd eigenfunctions. As noted beforehand, they vanish at the origin {x = 0) and are not affected by the 
presence of the Dirac delta source. Therefore, for the odd states the solution of the quantum harmonic oscillator is 
simply obtained. 

The problem of even states is more pronounced as the ground state is nodeless and even. A detailed investigation 
of the differential equation (A3) reveals that the even solutions are also written in terms of M and t/, but the 
initial conditions need to be chosen properly to account for the presence of the Dirac delta distribution. To find the 
required initial condition one integrates both sides of Eq. (A3) on a small interval around the origin, i.e., (—e,-l-e). 
Subsequently, all terms that cannot contribute in the small e limit are dropped. One can easily estimate which terms 
are significant by noting that both solutions of the homogeneous equation are regular around x = 0. Finally, the 
following result is obtained 

/'(+£)-/'(-£) = 25/(0). (A5) 

Let us recall that around x = 0 the solutions of the homogeneous equation behave for x > 0 like: 

U{-e/2,l/2,x^) = ^^ + Oix), (A6) 

U'i-e/2, 1/2, x^) = + 0(x), (A7) 

M(-e/2,l/2,x2) = l + 0(x2), (A8) 

M'(-e/2,l/2,x^) = -2ex-ke>(x^). (A9) 


The functions M are unable to satisfy the above initial condition. Therefore, we must pick up the solution expressed 
through the Tricomi functions U and from Eq. (A5) we obtain 


1 r(g^) 

~9 er(i^)- 


(AlO) 


The above expression cannot be solved with elementary methods. Nonetheless, for a given g the solution is straight¬ 
forward to obtain numerically. Finally, the exact wavefunction of the Hamiltonian (Al) is given by 


4'(x, X) = CnC, 17(-e/2,1/2, x"). 


(All) 


where Cn and Ce are the normalization constants for the x- and X-dependent portions, and the corresponding total 
energy of the system is simply A = e -|- n -I- 1 . 

Having the exact wavefunction of the system, we can analyze the possible cusp-like conditions at the particles 
coalescence points (x = 0). The wavefunction around x = 0 behaves as 

vl/(x,X)(x 1-2 ^1^^"^^^ \x\+0{x^), (A12) 

r(-2e) 

which is a direct consequence of the properties of U. The above expression can considerably be simplihed with the 
help of Eq. (AlO) which leads to 

'I'(x, X) oc 1 -I- 5 jxj -I- O(x^). (A-13) 

Strikingly, the behavior of the exact wavefunction around x = 0 depends only on the value of g. This is also the 
counterpart of the Kato’s electron-electron cusp condition. In analogy, Eq. (A13) can be written as 


d'i>ix,X) 


dx 


a;=0± 


±5VI/(0,X), 


(AM) 
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where the equality sign strictly holds for any X. Due to the requirement of even parity, the exact wavefunction 
possesses a derivative discontinuity at a; = 0 and a “cusp” at the origin. This discontinuity can be a considerable 
difficulty from the practical point of view. When the wavefunction is modeled with smooth basis set functions, the 
cusp may be very difficult to reproduce. 


Appendix B: Energy convergence with the basis set size - truncated exact wave function 

In the two-body case the expansion of the exact wavefunction in one-dimensional harmonic oscillator eigenfunctions 
can be written as 


^{x,X) = EE Cmn (pn{X). (Bl) 

m—0 n—0 

By recalling the form of the exact wave function, Eq. (3), one finds that the X-dependent part of the wavefunction 
is automatically described exactly. Therefore, the actual task in the two-body calculations is to approximate the 
x-dependent part of the wavefunction, by the solutions of the quantum harmonic oscillator eigenproblem, i.e., 

OO 

cl){x) = Ce e“5^V(-e/2, 1/2, =^'^Cn ifnix). (B2) 

n—0 

Due to the orthonormality of <fn{x) the coefficients c„ obey the formal relationship 

/ -l-oo 

dxe~i'^ U{-e/2,1/2,x'^)ipn{.x). (B3) 

-CX) 

Note that since the exact wavefunction is even, the coefficients vanish for odd n by symmetry. 

As the basis set, Eq. (2), is complete in the second Sobolev space, we can construct systematic approximations to 
the exact wavefunction by terminating the expansion (B2) at some n;,. This gives a family of approximants 


rib — l 

(j)M(^x) = ^ Cnipnix), 


n—0 


(B4) 


which are not normalized to the unity. The energy connected with for each Ub is given by 




(B5) 


Since the exact wavefunction, (/(x), is normalized, the exact energy of the relative motion is simply e = {(/IHxl/))- Let 
us also define the complementary function, 4'c^'’\x), given for each nt by the expression 


OO 

(j)l/"’\x) = <j){x) - ^ Cnipnix). (B6) 

n—Uh 


We can now precisely state that we are interested in finding the asymptotic expansion of en^ — e as —>■ oo. At this 
point we must stress that the energy is not strictly equal to the energy calculated in the same basis set. In fact, 
in the actual finite basis set calculations, the coefficients c„ are determined variationally rather than by projection 
onto the exact wavefunction. Therefore, the variationally determined coefficients have an opportunity to relax and 
accommodate to the incompleteness of the basis set. However, for large nj, this relaxation effect is small and virtually 
limited to the last coefficient, as shown by Carroll for the helium atom [117]. Therefore, the coefficients c„ and 

the energies are expected to be very close to the corresponding variational values. 

Let us recall Eq. (B3) and insert the explicit form of <fn(x) 


C, 7r-i/4 

Cn — , Xn, 

V2”n! 

f + °° 2 

Vn= dxe~^ U{—e/2,1/2, x^) Hn{x). 

J —OO 


(B7) 
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To evaluate the integrals we need to recall a particular integral representation of t/(a, b, z) 

Uia,b,z) = ^ / + 

r(a) Jo 


(B 8 ) 


which is valid for a > 0. Therefore, we confine ourselves here to the regime e < 0 which roughly corresponds to g < 0 
(attractive potential). However, one can show that the final result derived here remains valid also for e > 0. By 
inserting the above expression into the definition of v„ and exchanging the order of integrations, one finds 


v„ 


1 

rPP 



dtt 




0-\-OQ 




(B9) 


Let us consider only the inner integral for a moment, denoted by I . The change of variables u = x\/\ + t gives 

2 ^+oo 


I = 


\/l + ^ J —c 




VI T ^/ 


(BIO) 


To simplify the Hermite function under the integral sign one can make use of the following relation 


Hn{iu) 


nl2 





^Tfn- 2 *(X), 

7 I 


(Bll) 


where 7 = (1+t) ^. Upon inserting into the integral I and making use of the orthogonality of the Hermite polynomials 
one arrives at 


/ = V^7(7^ 




n! 


(B12) 


By returning to the definition of and slightly rearranging the following formula is obtained 

_ V^(-l)’"/2 n\ f°° 

r(-ie) {n/ 2 )\J, + 

and the remaining integral is elementary, so that 

_ V^(-l )”/2 n! 2 
r(-ie) (n/ 2 )!n-e' 

Therefore, the expression for c„ reads 

^ 7ri/4(-l)"/2 V^ 2 

' r(-ie) 2^/^{^ny.n-e' 


(B13) 


(B14) 


(B15) 


and the coefficients vanish for odd n, i.e., C 2 n+i = 0. The above expression is fairly difficult to handle because of the 
presence of the factorials. Fortunately, we require only their asymptotic forms, i.e. approximately valid for large n. 
With the help of the Stirling formula, n! —>■ VSttu one arrives at 


^2n 


. ^ 2 (- 1 )" 1 

'^r(—ie) 2n — e "nJ!^ ’ 


as n —>■ 00 . 


(B16) 


The first ingredient required for the presented derivation is the asymptotic formula for the square of the norm, 
for large value rif,. Taking advantage of the orthonormality of and normalization of the exact 
wavefunction one hnds 


= ^ cX (B17) 

n—n\, 

Since the second term of the above expression vanishes for large rib, one can write that —>• 1 as rib —>■ 00 

which is entirely sufficient for the present purposes. One can also verify that a somewhat more accurate formula, 
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accounting for the next term in the asymptotic expansion, would include a term proportional to However, 

this term does not contribute to the final results and is omitted for brevity. 

Because of the approximation derived for the denominator in Eq. (B5) one can rewrite it for large rib as 


By recalling the definition (B6) and rearranging one gets 

en. - e « . 

The first of the above matrix elements can be expanded to give 

f^cl(n+^]+g 

n—Tib ^ '' 

Similarly, for the second matrix element of Eq. (B19) one obtains 

OO / 1 \ oo 

d \n +-) + g <j){0) Y Cn^finiO), 

n—rii. ' ^ n—rii. 


n—rifj 


(B18) 

(B19) 

(B20) 

(B21) 


by noting that {(l>\<Pn) = Cn and ^(0) does not depend on rib- Note that the following two infinite sums are necessary 
for the evaluation of Eqs. (B21) and (B20) 


TN-J2^n(n + ^)= Y cL 2n + - , 


n—Ub /2 


= X! = X] C2n</?2n(0). 

n=ni, ni,/2 


(B22) 

(B23) 


The simplest way to evaluate these sums for large rib is to exchange the summation indices to run only over even 
values of n and subsequently insert the asymptotic formula for C 2 n, Eq. (B16). The resulting expressions are (showing 
only the first term of the Stirling series) 


T; 


( 1 ) 




N 


Ti-ur 


E 


T; 


( 2 ) 


2'"'^ n—nbl2 

2a 


N 


:i7. E 


1 2n+ ^ 

{2n — e)2 dn 

1 1 




where we have additionally used the relationship 


(— 1 )” 

</52n(0) —>■ —as n —>• oo. 


(B24) 

(B25) 

(B26) 


easily obtained from Eq. (2) by using the Stirling approximation. Finally, the above sums (taking into consideration 
all relevant terms in the Stirling series) can be evaluated with the help of the Euler-MacLaurin resummation formula 
which leads to 


( 1 ) ^ 2^/^C! / 1 , 86 + 7 

r(-ic)2 I 1/2+48 3/2 


O 


_ 

^ V^r(-ie) UE ' 48nE 


2 ' \"'b ’^'"6 , 

23/2(7, f 1 8e+ll' 


(» 7 ). 

o("7). 


Upon reinserting these expressions into Eqs. (B21), (B20), and (B19) one arrives at 

23/2(72 


- e = 


r(-ie)2 I y/n^ ' TT rib 


‘ +171) 


'(»r). 


(B27) 

(B28) 

(B29) 
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Let us also mention that with the same methodology as presented above one can derive further terms in the asymptotic 
expansion of — e. For example, the expression including additionally terms proportional to reads 


25 / 2(^2 I' I g^/2 

e™, - e - — 


15 + 8e\ 
48n^/V 


+ 0(n-2) . 


(B30) 


Note that the convergence proportional to Ij^/rTb in the leading order is extremely slow. The latter convergence 
formula may be employed to describe numerical results in the weak and intermediate coupling regime. In the next 
appendix we will instead outline a different approach, which yields a convergence formula valid for all coupling 
strengths. 


Appendix C: Energy convergence with the basis set size - variational wave function 

In Appendix B we have considered the problem of the convergence of the calculations with the increasing basis 
set size for the two-particle case by calculating the energy given by expectation value (B5) with the truncated exact 
wave function (B4) for a hnite interaction strength g. Now, we will present an alternative derivation, allowing for 
a variational relaxation of the coefficients Cn in the truncated wave function (B4). Obviously, since the functional 
that we minimize is quadratic, the resulting wave function corresponds to the exact ground state wave funetion of the 
truncated (projected) Hamiltonian. In the following, we shall first consider the numerical convergence of the inverse 
coupling constant versus Ub, at fixed energy e, and finally derive the convergence of the energy at fixed coupling. 

By minimizing the expectation value of the Hamiltonian on the wavefunction (B4) with respect to its coefficients 
Cn, one finds 


0 = ^ {<f\Ha; - e\ip) = (n - en„)Cn + gfn^Cn' fn 


(Cl) 


where 2i/^/„ = = ‘Pn(O) is the n-th harmonic oscillator wavefunction for the relative coordinate, evaluated 

at X = 0. We have / 2 n+i = 0, and 


1 (2n-!)!!_ 1 (2n)! 

2n!! “ v^22-(n!)2 ’ 


(C2) 


By solving Eq. (Cl) for g in a procedure similar to that of Busch et al. [78], one finds the result presented in Eq. (4) 
of the main text. 


y2er[(l - e)/2] ■ 


(C3) 


Let us now define x(e) = g ^(e). In numerical calculations, one has to truncate the basis to a certain maximum 
number of states nb, and therefore obtains only an approximate value, 


Xn.ie) 


1 

5n„(e) 


rib —1 

E 


n—O 


(/2n)^ 

e — 2n 


(C4) 


The convergence rate of this formula at a fixed energy e, as a function of the number of states included in the 
summation, Ub, is given by 


f \ / \ (/2n)^ 

*„(e)-x(e)= 

71—Ub 


1 / 1 5 -b 4e \ 


(C5) 


as may be found using Stirling’s formula, and then expanding the fraction inside Eq. (C3) for 2n ^ jej. 

Let us now consider the numerical error Ae„^ = — e obtained by computing the energy at a fixed g with a basis 

set containing up to and including Ub functions. The error is given by the solution of the implicit equation 


=a;„,(e„J, 


(C6) 
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which can be rewritten as 


Aa;„^ (e + ) = -[a;(e + ) - x{e)\ . 


(C7) 


On the left-hand-side of the above equation we can now use Eq. (C5). Assuming that Ae„j can be expanded in powers 
of Ijy/ub and equating the coefficients on both sides, we obtain 


A- 


\ 1 

Ax” 

1 

5 + 4e 


fx”V x'” 

1 1 


2x' 

Ub 
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+ 2 

\x') 3a;' 

-f i 


+ 0[n-^), 


(C8) 


with A = —1/{'\/2'kx'), and x' = df\g ^(e)]. In the vicinity of the Tonks-Girardeau point where e ~ 1 we have 
e = 1 - so that A = 2/7r^/^. In the weak coupling limit e ~ 0 instead, we have e = gj^f^ -I- ... so 

that A = g^ j Both these limits coincide with the analytical results of Ref. [53]. 
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